C-------------------------------------------------------------------------------
C  x: Ax = y
C-------------------------------------------------------------------------------
      subroutine cg(prec, Nx, Ny, Nz, A, y, tol, x, xreal)
      implicit none
      
      integer Nx, Ny, Nz, prec
      
      double precision A(Nx,Ny,Nz,4), y(Nx,Ny,Nz), tol
      double precision x(Nx,Ny,Nz),  xreal(Nx,Ny,Nz)
      
      double precision, dimension(:,:,:), allocatable :: p, b, r, r_old,
     1                                                   D, temp_vec, z
      double precision alpha, temp
      double precision tol2
      double precision rz, rz_old
      double precision norm_res2, norm_res2_0
      
      integer Nt, i, MAXITER
      
      parameter(MAXITER = 1000)
      
      
      allocate(p(Nx,Ny,Nz))
      allocate(b(Nx,Ny,Nz))
      allocate(r(Nx,Ny,Nz))
      allocate(r_old(Nx,Ny,Nz))
      allocate(D(Nx,Ny,Nz))
      allocate(temp_vec(Nx,Ny,Nz))
      allocate(z(Nx,Ny,Nz))
      
      Nt = Nx*Ny*Nz
      
      tol2 = tol*tol
      
      call make_cholinc(Nx, Ny, Nz, A, D)
      
      call fill_vec(Nt, 0d0, x)
      call copy_vec(Nt, y, r)
      
      call precond(prec, Nx, Ny, Nz, A, D, r, z)
      
      call copy_vec(Nt, z, p)
      
      call dot_prod(Nt, r, r, norm_res2_0)
      call dot_prod(Nt, r, z, rz)
      
      do i=1,MAXITER
	
	call matvec_prod(Nx, Ny, Nz, A, p, temp_vec)
	call dot_prod(Nt, p, temp_vec, temp)
	alpha = rz/temp
	
	call lin_comb(Nt, x, p, alpha, x)
	call lin_comb(Nt, r, temp_vec, -alpha, r)
	
	call dot_prod(Nt, r, r, norm_res2)
	norm_res2 = norm_res2 / norm_res2_0
	
	write(6,'(I10,2X,E9.3)') i, sqrt(norm_res2)
	
	if (norm_res2 .le. tol2) goto 100
	
        call precond(prec, Nx, Ny, Nz, A, D, r, z)
	
	rz_old = rz
	call dot_prod(Nt, r, z, rz)
	
	call lin_comb(Nt, z, p, rz/rz_old, p)
      end do
      
      write(6,"('CG didn''t converge')")
      
  100 continue
  
      call lin_comb(Nt, x, xreal, -1d0, temp_vec)
      call dot_prod(Nt, temp_vec, temp_vec, temp)
      temp = temp / norm_res2_0
      write(6,"('Exact error:')")
      write(6,'(I10,2X,E9.3)') i, sqrt(temp)
  
      deallocate(p)
      deallocate(b)
      deallocate(r)
      deallocate(r_old)
      deallocate(D)
      deallocate(temp_vec)
      deallocate(z)
  
      return
      end
C-------------------------------------------------------------------------------
C  
C-------------------------------------------------------------------------------
